home *** CD-ROM | disk | FTP | other *** search
/ FM Towns: Free Software Collection 9 / FM Towns Free Software Collection 9.iso / t_os / tool / pcom / src / fft.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-16  |  1.7 KB  |  99 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include <nslib.h>
  4. #include <pc_inc.h>
  5.  
  6. static int tw[N], s[2][N], c[N], ex1[N/2], ex2[N/2];
  7.  
  8.  
  9. void sin_cos( void )
  10. {
  11.     double p=0, d=3.1415927*2/N;
  12.     int i;
  13.  
  14.     for( i=0; i<N; i++ )
  15.         {
  16.         s[0][i] = (int)(sin( p )*65536);
  17.         s[1][i] = -s[0][i];
  18.         c[i] = (int)(cos( p )*65536);
  19.         tw[i] = 35389 - ( (30146*c[i] +0x8000 ) >> 16 );
  20.         p += d;
  21.         }
  22. }
  23.  
  24.  
  25. void bitrev( void )
  26. {
  27.     int i, j, ii, bit, di=0;
  28.  
  29.     for( i=1; i<N-1; i++ )    /* bit reversal */
  30.         {
  31.         for( j=0, ii=i, bit=0; j<EXN; j++, ii>>=1 )
  32.             {
  33.             bit <<= 1;
  34.             bit |= (ii&1);
  35.             }
  36.         if ( i < bit )
  37.             {
  38.             ex1[di] = i;
  39.             ex2[di++] = bit;
  40.             }
  41.         }
  42.     ex1[di] = 0;
  43. }
  44.  
  45.  
  46. void time_window( short *buf )
  47. {
  48.     int i;
  49.  
  50.     for( i=0; i<N; i++ )
  51.         buf[i] = (short)( ( (int)buf[i] * tw[i] +0x8000 ) >> 16 );
  52. }
  53.  
  54.  
  55. /* inv = 0 ＦＦＴ    inv = 1 逆ＦＦＴ */
  56. void fft( short *rl, short *im, int inv )
  57. {
  58.     int i, j, k, j1, j2, tim, lenb=N, numb=1, w;
  59.     short xr, xi, yr, yi;
  60.  
  61.     for( i=0; i<EXN; i++ )
  62.         {
  63.         tim = 0; lenb /= 2;
  64.         for( j=0; j<numb; j++ )
  65.             {
  66.             w = 0;
  67.             for( k=0; k<lenb; k++ )
  68.                 {
  69.                 j1 = tim+k;  j2 = j1+lenb;
  70.                 xr = rl[j1]; xi = im[j1];
  71.                 yr = rl[j2]; yi = im[j2];
  72.                 rl[j1] = xr + yr;
  73.                 im[j1] = xi + yi;
  74.                 xr = xr - yr;
  75.                 xi = xi - yi;
  76.                 rl[j2] = (short)( ( (int)xr*c[w] - (int)xi*s[inv][w] +0x8000 )
  77.                         >> 16 );
  78.                 im[j2] = (short)( ( (int)xr*s[inv][w] + (int)xi*c[w] +0x8000 )
  79.                         >> 16 );
  80.                 w += numb;
  81.                 }
  82.             tim += (lenb*2);
  83.             }
  84.         numb *= 2;
  85.         }
  86.  
  87.     for( i=0; ex1[i]; i++ )    /* bit reversal */
  88.         {
  89.         k = rl[ex1[i]]; rl[ex1[i]] = rl[ex2[i]]; rl[ex2[i]] = k;
  90.         k = im[ex1[i]]; im[ex1[i]] = im[ex2[i]]; im[ex2[i]] = k;
  91.         }
  92.  
  93.     for( i=0; i<N; i++ )
  94.         {
  95.         rl[i] /= SQRTN;
  96.         im[i] /= SQRTN;
  97.         }
  98. }
  99.